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ABSTRACT 

Realistic numerical simulations, i.e., those that make minimal use of ad hoc 
modeling, are essential for understanding the complex turbulent dynamics of 
the interiors and atmospheres of the Sun and other stars and the basic mech¬ 
anisms of their magnetic activity and variability. The goal of this paper is to 
present a detailed description and test results of a compressible radiative MHD 
code, ‘StellarBox’, specifically developed for simulating the convection zones, 
surface, and atmospheres of the Sun and moderate-mass stars. The code solves 
the three-dimensional, fully coupled compressible MHD equations using a fourth- 
order Pade spatial differentiation scheme and a fourth-order Runge-Kutta scheme 
for time integration. The radiative transfer equation is solved using the Feautrier 
method for bi-directional ray tracing and an opacity-binning technique. A spe¬ 
cific feature of the code is the implementation of subgrid-scale MHD turbulence 
models. The data structures are automatically configured, depending on the 
computational grid and the number of available processors, to achieve good load 
balancing. We present test results and illustrate the code’s capabilities for sim¬ 
ulating the granular convection on the Sun and a set of main-sequence stars. 

The results reveal substantial changes in the near-surface turbulent convection 
in these stars, which in turn affect properties of the surface magnetic fields. 

For example, in the solar case initially uniform vertical magnetic fields tend to 
self-organize into compact (pore-like) magnetic structures, while in more mas¬ 
sive stars such structures are not formed and the magnetic field is distributed 
more-or-less uniformly in the intergranular lanes. 
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Subject headings: magnetohydrodynamics (MHD) - radiative transfer - radia¬ 
tion: dynamics - convection - turbulence - methods: numerical - Sun: atmo¬ 
sphere, interior, magnetic fields - stars: atmospheres, interiors, magnetic held 


1. Introduction 


During the last few decades, there has been a major increase in the use of high- 
performance computing in computational physics in general and in space science in par¬ 
ticular. Massively parallel algorithms have been developed to implement accurate physical 
models in order to provide efficient and realistic simulations of astrophysical phenomena, for 
instance, stellar and solar interior dynamics, which is the focus here. Such computational 
tools are used to support and analyze ground and spacecraft observations. The coupling 
between observations and numerical simulations is necessary for improving our understand¬ 
ing of complex phenomena on the Sun and other stars because the two approaches are 
complementary. Recent high-resolution observations of the Sun have uncovered a rich small- 
scale dynamics which plays a key role in the observed large-scale phenomena. For example, 
observations and simulations led to the discovery of small-scale vortex tubes generated in 


intergranular lanes ( 

Brandt et al. 

1988; 

Bonet et al. 

2008, 

2010; 

Steiner et al. 

2010 

), which 

are a source of the Sun’s acoustic emission and magnetic flux concentrations ( 

Kitiashvili 


et al. 2010, 2011, 2012). Stellar observations from the Kepler mission and ground-based 


telescopes have revealed that solar-type acoustic oscillations are common in other stars, giv¬ 
ing rise to the rapid development of asteroseismology. Also, these observations have found 
a broad range of magnetic activity, including magnetic activity cycles, starspots, and sur¬ 
prisingly strong energy release events (‘superflares’). Because the structure and dynamics of 
the surface of stars other than the Sun are not resolved in observations, it is important to 
develop numerical simulations capable of reproducing realistic stellar conditions. For solar 
conditions such simulations can be verified by comparing with high-resolution observations, 
providing confidence for using these simulations in the interpretation of stellar observations. 

One particular class of these numerical investigations is the realistic simulation of the 


upper part of the convection zone and the lower part of the atmosphere of the Sun (Nordlund 

1982; Keller et al. 

2004; Garlsson et al. 2004; Vogler et al. 2005; Stein & Nordlund 2006; 

Georgobiani et al. 2012) and other stars (Giorgobiani 2003; Steffen et al. 2005; Kupka 2009; 

Steffen et al. 2009; 

Muthsam et al. 2011J Georgobiani et al. 2012; Ludwig & Kucinskas 2012; 

Grimm-Strele et al 

. 2015; Mundprecht et al. 2015). Such investigations typically consist of 


compressible radiative MHD simulations in which the governing conservation equations of 
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mass, momentum, energy, and magnetic flux are integrated in time. The stellar composition 
is modeled using a specific tabular equation of state (EOS) and optical opacity, both based 
on given stellar chemical abundances. The computational cost of this approach is high, and 
currently only a relatively small part of a stellar body can be simulated with a resolution 
sufficient to study the turbulent dynamics of the surface and atmosphere in detail, e.g., to 
resolve in detail granulation and acoustic events. Nevertheless, these small-scale simulations 
provide knowledge about the structure and dynamics of the near-surface turbulent layer— 
one of the most complex regions inside stars. This layer plays a key role in energy and mass 
transport from a star’s interior to its corona and wind, as result of strong coupling among 
plasma dynamics, radiation, and magnetic fields. This leads to turbulent convection form¬ 
ing high-speed downdrafts driven by radiative cooling and to the magnetic field becoming 
organized into strong-field structures that dominate the outer envelope. Knowledge of this 
dynamics, gleaned from high-resolution local simulations, can then be used in global-star 
models with relatively low resolution. 


The current paper describes the implementation and testing of the radiative MHD code 
‘StellarBox’. The code provides realistic simulation of solar and stellar convection zones and 
atmospheres and has been used for a wide range of problems, such as multi-scale dynamics 
and self-organization processes in turbulent magneto-convection, acoustic wave excitation, 
formation of stable magnetic structures (such as pores and sunspots), eruptions, generation 
of magnetic fields by local dynamos, simulation of specific local conditions (e.g., sunspot 
umbrae and penumbrae), and interactions of the turbulent surface and subsurface with 
atmospheric layers (Jacoutot et al. 2008; Kitiashvili et al. 2009, 2010, 2013a) . In addition, 
the code makes it possible to use as initial conditions various models of interior structure 
pre-calculated for stars with specific chemical compositions, masses, and rotation rates. For 
F- and A-type stars it is feasible to extend the computational domain into the radiative zone 
to study the dynamics of convective overshoot layers between the radiative and convection 
zones. In this paper we first focus on a description of the code and tests of its numerical 
methods. Then, as an example application, we show the code’s capabilities in simulating 
turbulent magnetoconvection on the Sun and several main-sequence, moderate-mass stars. 
A detailed analysis of the stellar simulations will be presented in a separate paper. 


The paper is structured as follows. After a brief description of the code modules and 
components, the governing equations and the physical models are given. Details of the 
numerical methods and boundary conditions are described. After this, validation and seal- 
ability tests are provided. Finally, example results from simulations of solar and stellar 
magnetoconvection are presented. 
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2. Code Description 


StellarBox was developed at the NASA Ames Research Center to solve the three- 
dimensional, fully coupled compressible radiative MHD equations using a fourth-order com¬ 
pact implicit Pade scheme for spatial differentiation and a fourth-order Runge-Kutta scheme 
for time integration. The choice of the numerical scheme is motivated by previous experi¬ 
ence in numerical turbulence modeling (e.g., Miyauchi et al.||1993 ; Chow fe Moin||2003 ). The 
computational domain is a rectangular volume typically encompassing a horizontal span of 
tens to hundreds of megameters and a vertical span of a few to several tens of megameters. 
The domain is discretized using a Cartesian grid. Spatial resolution is typically in the range 
of 100 km down to 6 km in the horizontal directions. Arbitrary (user-specifiable) mesh 
stretching is used in the vertical direction. The vertical grid spacing typically increases with 
depth and ranges from from 10 km to 100 km. Periodic boundary conditions are used in 
the horizontal directions and characteristic boundary conditions are applied in the vertical. 
Stellar rotation is accommodated by an f-plane approximation at a user-specified latitude 
and rotation rate. 

The code incorporates the following user-selectable subgrid-scale turbulence models for 


the transport of heat and momentum: compressible Smagorinsky (1963) model, dynamic 


Smagorinsky ( Moin et al.||1991 ; Germano et al. 1991), and implicit hyperviscosity (Caughey 
& Jothiprasad 2002). 


The radiative transport equation (RTE) is solved at every full time-step by using a ray¬ 


tracing, long characteristic algorithm along 18 rays. A second-order Feautrier (1964) method 


is used to solve the RTE along each ray. Frequencies are opacity-binned into logarithmic 
bins, typically four in number. At depths at which the medium is optically thick, diffusive 
radiative transport may be optionally used instead of ray-tracing for improved efficiency and 


accuracy. A tabular Equation of State (EOS) (Rogers et al. 1996) and opacity are used. 


StellarBox is a massively parallel code that uses algorithms optimized for parallelization. 
The parallel data structures are automatically configured, for given mesh dimensions and 
number of processors, to optimize load balancing. Differentiation and the radiative transfer 
solution are accomplished by transposing these data structures so that x, y, or z is memory- 
resident as needed. 

The overall code structure consists of four main components: Radiation, Time Advance, 
Utility, and Thermodynamics. Figure [l] shows the code components and subcomponents. 
Their functions are described in the following bullets: 


• Radiation: implements the opacity-binning model, interpolation algorithms for the 
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opacity and radiation source functions, and the RTE solver. 

• Time Advance: implements the convective and diffusive flux calculations, the shock¬ 
capturing and turbulence models, and the boundary conditions to perform the time 
advance. 

• Utility: implements various utility functions needed by the other modules. 

• Thermodynamics: implements the equation of state (EOS) and calculates thermo¬ 
dynamic properties. 


2.1. Governing Equations 

The equations governing compressible radiative MHD flows are conservation of mass, 
momentum, energy, and magnetic flux. The code solves grid-cell averaged conservation 
equations for these quantities. 

The mass conservation equation: 

dtp + (puj)j = 0 ( 1 ) 

where d t denotes the time derivative operator. The subscript , j denotes the space derivative 
operator in the j th direction with j = {1,2,3}. Quantity p is the averaged mass density 
and Uj is the Favre-averaged, i.e., density-weighted average, velocity component in the j th 
direction. 

The momentum conservation equation with gravitational and magnetic fields can be 
written, in the stellar rotation frame, as: 

df (pUi) T (pUiUj T pSij) j ffq'.j T fdij.j P0.i 26,^/. ffljpu k (2) 

where p is the pressure, 8ij the Kronecker symbol, (ft the gravitational potential, n, 7 the 
viscous stress tensor, the permutation tensor, Q :] the stellar mean-rotation vector, and 
Bij the magnetic stress tensor 

Bij = —BiBj - —B k B k 8ij (3) 


The total energy conservation equation reads 


d t E +[(E + p)u j ] j = 


j T (11 ij Ui ) j T (B i: jii i ) j 


( 4 ) 
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where II is the viscous tensor, Q is the non-radiative heat flux (diffusive and Joule heating), 
Q T j ad is the radiative heat flux, and E is the total energy per unit volume given by: 

E = pe + —piijUj + p(f> + ——BjBj (5) 

2 o7r 

where e is the internal energy per unit mass. The radiative heat flux is obtained by solving 
the radiative transfer equation, and the viscous tensor and the non-radiative heat flux are 
written in terms of transport coefficients and gradients of resoived variables. 


The magnetic flux conservation equation can be expressed in Gaussian units as 


fit Bj T ( Uj Bi UiBj) ■ 


47ra 


\B-i,j B 




( 6 ) 


where Bj is the magnetic flux density in the j th direction, c the speed of light, and <r the 
electrical conductivity. 


2.2. Subgrid Stress Tensor 


Since it is currently impossible to achieve realistic Reynolds numbers in solar and stellar 
simulations, StellarBox utilizes a Large-Eddy Simulation (LES) approach in which models 
are used to approximate the effect of subgrid motions on the resolved scales. Modeled subgrid 
quantities will be denoted by a subscript T. The subgrid stress tensor is represented in the 
form: 


fl ij -pT ( Sij ^ Llk kSjj ) BfSjj 


( 7 ) 


in which px is the subgrid viscosity given by the Smagorinsky (1963) model, augmented for 
compressible flows with a shock capturing term: 

p T = pA 2 (Cs\S\ + Cc\uk t k\) (8) 

where Cs is a Smagorinsky coefficient, Co is a shock-capturing coefficient, A is the grid 
spacing, (S'! = \J{2SijSij), and is the strain-rate tensor 

Sij = — (Ui,j + u j,i) j (9) 

kr is the subgrid kinetic energy density given by 

k T = ^C c pA 2 \S\ 2 (10) 

where C'c is a second Smagorinsky coefficient. The code provides a dynamic option for deter¬ 


mining Cs and Cc based on the associated Germano et al. (1991) identities. Alternatively, 
they can be specified by the user, as is Co- 
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2.3. Heat Fluxes and Turbulent Resistivity 

The total non-radiative heat flux is the sum of two contributions: a term due to tem¬ 
perature gradients, that is, Fourier’s law, and a second due to Joule heating: 

Qj = KtT,) + ( 4 ^) — BjBjj ), ( 11 ) 

where kt is the subgrid heat conductivity 


c p fl T 

KT ~1W 

in which Prx is the turbulent Prandtl number and c p is the specific heat at constant pressure 
per unit mass taken from the Equation Of State (EOS) tables. The turbulent Prandtl number 
Prx can be specified by the user; it is typically taken to be near unity. 



The turbulent electrical conductivity, ay, in Gaussian units is given by 

c 2 

(Jj 1 


Aitt/t ’ 


(13) 


where c is the speed of light and r/x is the turbulent magnetic diffusivity, modeled following 
Balarac et al. (2010) and Theobald et al. ( |1994| ): 


r/T — C#A 


2 \ e ijkBk,j\ 


\fp 


(14) 


where Cb a user-specified constant chosen through comparison to fully resolved MHD tur¬ 
bulent simulations. A typical value is 0.25. 


2.4. Radiative Heat Flux 

The radiative heat flux is given by 




f lj I u (Cl, x) (fndu, 


(15) 


where the quantity H is a direction unit vector, u is the frequency, and I u is the radiative 
intensity at frequency u given by the radiative transfer equation: 
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n ■ vj „ (to, x, t) = Xu (x, t) [s u (x, t) - i v (to, x, t)] , 


(16) 


where non-isotropic scattering, polarization, and the finiteness of the speed of light are 
neglected. S v is the radiation source function and Xu the opacity. Eq. (16) is solved using 
the Feautrier (1964) method, which is based on bi-directional ray tracing. The algorithm 
implemented in the code uses 18 rays (9 bi-directional rays), shown in Fig. [2j 


In the Feautrier method, the intensity function along a ray is divided into forward- and 
backward-travelling intensities (/+, I~ respectively) as a function of distance s along the 
ray: 


dl + 

- = Xu (s v - /+) 


ds 

dl 

ds 


Xu I u ) 


By subtracting Eq. (18) from Eq. © , we obtain 

ax-i, 


ds 


= 2 XuS v - Xu (It + It) 


and by adding Eq. © and Eq. ( fl8| ) , we have 

d (/+ + /“) 


ds 


= -x , X - K) 


(17) 

(18) 


(19) 


( 20 ) 


Defining If = It ~ I v and /^ = /+ + I v and combining Eqs. 
single second-order differential equation 


(19) and (20), we obtain a 


(Plt_ _ 1 dxv dlt 

ds 2 Xu ds ds 


- xX = -ixlSv 


( 21 ) 


Equation (21), as written, is for a monochromatic intensity / ; f. As mentioned earlier, 
the code treats frequency space by using an opacity-binning technique, in which the opacity 
of a given bin b, Xb, is taken to be the Rosseland mean opacity of the frequencies assigned 
to that bin. Similarly, the source function Sb for bin b is the source S u integrated over those 
frequencies. The RTE for the integrated intensity A in bin b is then 


d 2 /f i dxbdijg 
ds 2 Xb ds ds 


Xblb ~ — 2 XbSb 


( 22 ) 
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Eq. (15) becomes, expressed in terms of bins 


Tib p 

Q- ad = / n, / 6 (o, x ) d 2 Q 

6=1 ^ 47r 


(23) 


where rib is the number of opacity bins, typically four in running StellarBox. Upon discretiz¬ 
ing Q. the integral over Q in Eq. (23) becomes an appropriately weighted sum over the 18 
rays shown in Fig. [2j The weights are chosen to maximize the Taylor-series order of accuracy 
of the solid-angle integration. 


2.5. Numerical Methods 

The computational domain is a rectangular volume of dimensions L x x L y x L z defined 
by the user. The volume is discretized into a cartesian grid, where the number of grid points 
along each axis, n x , n y , and n z , is set by the user as well. By convention, negative values of 
^ correspond to locations below the nominal photosphere and positive values above it. 

The basic finite difference differentiation scheme is 4th-order Pade and is used to com¬ 
pute derivatives appearing in the diffusive fluxes and also the derivatives of the convective 
and diffusive fluxes themselves in equations and (j6j). The scheme must of course 

be supplemented with boundary conditions, as described in the next section. 

The basic stencil to compute the first derivative T' of a quantity T in the j th direction 
U = {!, 2 ,3}) is 


. J h 


k -1 




\ T ' = 
4 k+1 


4 hi 


(J r k +1 - Fk-l ), 


(24) 


• th direction and h], is a measure of the 


where the subscripts are grid point indices in the j 
grid spacing at point k in that direction. Eq. (24) is all that is required in the periodic x 
and y directions, where 


h\ = Ax = — and hi = Ay = —. 


klnr 


n qi 


(25) 


In the vertical direction z (j = 3), the grid spacing is arbitrary, meaning that hf. may 
depend on k. Typically, the 2 -grid is stretched below the photosphere where structures are 


larger. For one-sided derivatives at the boundaries, Eq. (24) is supplemented, e.g., at the 
bottom in z, by the second-order boundary form 
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5(^1+^) = 4 (^ 2 -^) ( 26 ) 

and analogously at n z . Some ^-derivatives are computed by specifying derivative values at 
the boundaries, such as where characteristic boundary conditions are used. h\ is defined as 


hfc 2 (^fe+i %k— 1 ) V h \2,...n z 1} 

h\ = z 2 - zi 

l n z = Z n z ~ Z n z -1 


(27) 

(28) 
(29) 


The source term in Eq. (jd]) involves the computation of V • Q rad which, as mentioned 


in Sec. (2.4), requires solving Eq. (22). The latter is discretized using second order finite 


differences and solved for the quantity — 2S& so that the limiting case of an optically thick 
medium is handled accurately. 


The system of equations Eqs. (24) and the discretized form of Eq. (22) are tridiago¬ 


nal and are solved using the Thomas algorithm. The code has various tridiagonal solvers 
implemented for each type of boundary condition. 


2.6. Boundary Conditions 


In order to close the system of equations 0.0.0 , and g, boundary conditions are 
required. Periodic boundary conditions are applied in the horizontal directions (j = 1,2). 
For the z direction, the user has a choice of both boundaries closed (impenetrable, i.e., 
u 3 = 0), or both open, or the top open and the bottom closed. The open boundaries are 
implemented by a characteristic method (e.g. Sun et al. 1995). This top-open, bottom-closed 
set of boundary conditions is the one typically used in StellarBox. 


To simulate the energy flowing from the stellar core, the bottom boundary condition for 
the total energy is modified by adding an incoming energy flux per unit area equal to the 
stellar value. In addition, to preserve exact conservation of mass, momentum, energy, and 
magnetic flux in the domain, inward fluxes equal to the sum of those convected and diffused 
outward through the top and bottom boundaries are introduced at the bottom boundary as 
an areal average. This does not apply to the outgoing radiative flux at the top boundary — 
the system is allowed to find its own radiative equilibrium wherein the radiative flux emitted 
through the top boundary statistically balances the incoming energy flux at the bottom 







11 


boundary. The state that attains this last condition serves as a sanity check on the general 
correctness of the simulation. 


2.7. Time Integration 

The discretized system of equations ®, 0, and ([6j) can be written compactly as 


dXJ 

dt 


77(U) 


(30) 


where 77(U) includes all the spatially discretized terms, and U = (p, pu, E, B) T is the vector 
of conserved variables. Eq. (30) is solved using the following 4 th -ovdev Runge-Kutta scheme: 


U n+1 — U n H—— ( k\ + 2k 2 + 2k 3 + Aq) (31) 

b 

At is the time step, the superscript n stands for the time step, and the quantities Aq, k 2 , k 3 , 
and k 4 are defined as follows 


ki 

= 77(U n ) 

(32) 

k 2 

= K ( u " + fc 'T) 

(33) 

h 

= K ( u " + fe y) 

(34) 


= n(U n + k 3 At) 

(35) 


3. Parallel Scaling 

It is essential, for efficient parallel computation, to understand the scaling properties of 
a massively parallel code such as StellarBox. Various means of describing these properties 
have been devised. Here we present the results of our scaling tests in a particularly simple 
format based on comparison to an ideally performing code and computer system. By “ideal” 
we mean that the time associated with communication among the processors increases no 
faster than the amount of data communicated. Since all the significant computation in 
StellarBox involves an equal amount of work for each processor and for each grid point, an 
expression for the time, fideab for such an ideal system to compute one time step would have 
the following form: 
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^ideal ^ 


Tl x TlyTlz 


proc 


(36) 


in which n x ,n y , n z are the mesh dimensions and N proc is the number of processors. The factor 
a , which is the time per step per processor per grid point, is constant in the ideal case for any 
mesh dimensions or number of processors. Of course, any real code and computer system 
will show performance degradation, reflected by an increasing a, as the number of processors 
increases for a fixed problem size, principally due to contention for limited resources such 
as inter-processor communication hardware. In Figure [3] we show the a values obtained for 
various numbers of processors on two different mesh sizes: 1024 3 and 1500 3 ; all runs were 
conducted on the Pleiades computer system at NASA Ames Research Center. It is clear from 
this figure, for meshes of these approximate sizes, that the code behaves well — only slowly 
growing a — up to 45,000 processors or so, and that performance is seriously degraded for 
65,000 processors with a mesh of 1024 3 , on this computer system. Larger mesh dimensions 
would presumably result in better performance for 65,000 processors and beyond, though we 
have not tested this at this time. 


4. Test Cases 

SolarBox has been exercised for a large number of test cases to ensure physical and 
numerical accuracy. A subset of the more interesting tests is given in the subsections to 
follow. 


4.1. Sod Shock Tube 


The Sod shock tube problem (Sod 1978) is often used as a one-dimensional test case to 
check the ability of a compressible code to capture shocks, contact surfaces, and rarefaction 
waves present in the flow. The initial conditions are: 


, , f (0.125,0,0.1) if x <0.5 

(l.,o,l.) if £ > 0.5 

The gas is taken to be perfect with a specific heat ratio 7 = 5/3. Different runs were 
performed using different grid resolutions and different numerical diffusion coefficients in 
order to find a balance between damping numerical instabilities and resolving the sharp dis¬ 
continuities. In the case of the Sod shock tube, and as a general rule, less numerical diffusion 




13 


is required for ID, non-radiative, non-stratified problems than for stellar simulations. Figure 
[4] shows the numerical solution for pressure and density profiles along the x axis at time 
t = 0.1, using 480 grid points. The solid line is the numerical solution obtained using the 
differentiation, shock capturing, and time-advance algorithms in the code, and the symbols 


show the solution provided by an exact Riemann solver (Toro 1997). Good agreement is 


obtained for both the positions and the amplitudes of all the flow structures (shock, contact 
surface, and expansion fan). 


4.2. Orszag-Tang Problem 


The Orszag & Tang (1979) problem is a popular MHD test for two dimensional codes 


and is also known as the V • B = 0 test condition. It is used to check the robustness of 
a code in handling MHD shocks, including shock-shock interactions. The initial conditions 
are: 


u = — 
B x = 


25 

P ~ 36tt’ 

5 

P= l2n 


(37) 

-sin (27 ry); 

v = sin(27nr) 

Vx,ye [0,1] x [0,1] 

(38) 

-sin(2y ry); 

By = sin(47rx) 

\/x,y £ [0,1] x [0,1] 

(39) 

was performed using a 512 x 

512 grid. At time t = 0.25, 

Figure [5] 


shows snapshots of the density, magnetic energy, and kinetic energy fields. The images are 


directly comparable with ones in Stone et al. (2008) and good quantitative agreement is 


obtained even though the numerical methods used in the two codes are quite different. 

Figure [6] shows the pressure and density profiles for t = 0.25 at y = 0.4277. These 


results compare very well with those in Londrillo & Del Zanna (2000), Jiang & Wu (1999), 


and Ryu et al. (1998). 


4.3. Brio and Wu Shock Tube 


This test is a simulation of an MHD shocktube (Brio & Wu 1988); the hydrodynamic 


portion of the initial conditions are the same as for the Sod shock tube problem. However, 
the B field makes the algebraic equations of the Riemann problem highly non-linear and 
complex in a five-dimensional parameter space. Moreover, the presence of so-called non¬ 
regular waves in the MHD system causes the Riemann problem to be non-unique in some 
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cases (Torrilhon 2003). The right and left states are initialized as follows: 


(p,u,v,B y ,B z ,p ) t= o = 


(0.125,0,0,1,0,0.1) if x <0.5 
(1,0,0,-1,0,1) if x > 0.5 

where B x — 0.75 is constant and 7 = 2 . The numerical results are compared to the solution 
provided by the exact MHD Riemann solver (Exact-RS) of Torrilhon (2003). Figure [T] shows 
the density, pressure, and B y profiles at t = 0 . 1 . Good agreement is obtained for both the 
regular and non-regular waves. 


4.4. Radiative Transfer Test 


In order to test the ray tracing algorithm that is used in StellarBox to solve the radiative 
transfer equation (Eq. 16), two three-dimensional radiative transfer computations were 


performed using a box of 6 x 6 x 6 Mm that includes 1 Mm of the lower solar atmosphere; 
the grid dimensions were 64x64x128. In the first computation, the ray tracing algorithm 
was used everywhere; in the second, it was only used in the part of the domain defined 
by z > — 2 Mm (as previously mentioned, by convention negative values of z correspond to 
locations below the nominal photosphere and positive values above it), and a purely diffusive 
treatment was used everywhere else (z < — 2 Mm). The latter is a good approximation in 
optically thick regions, which is the case below z ~ — 2 Mm. The ray-tracing portion was 
computed using four opacity bins as usual, as described above, and the diffusive portion 
was computed using a single, all-frequency Rosseland-average opacity value at each point to 
obtain the diffusion coefficient. 


Figures [ 8 ] and [9] compare the local radiative cooling obtained for these two simulations. 
The profiles are at a single time and are averaged over x-y planes. They clearly agree very 
well over the full domain in z (Fig. [ 8 ]), and, in the optically thick region below z = —2 Mm 
(Fig. §, agreement is also good, though it is worth noting that the profile is smoother for 
the diffusive method, most likely because it does not involve a discrete angular quadrature; 
in any case, the values are very close between the two approaches. 


5. Simulations of Stellar Magnetoconvection 

In this Section we present some results of our simulations of solar and stellar magneto¬ 
convection. While the main goal is to demonstrate the code’s capabilities for understanding 
the dynamics of near-surface turbulent convection, the results give important insights into 
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stellar surface dynamics and magnetic effects. 


5.1. Effects of Numerical Grid Resolution 


Despite substantial growth in computing power, our computational capabilities for 3D 
simulations are still quite limited, and it is important to investigate the effects of numerical 
grid spacing on the resolution of the important turbulent convective structures, such as 


granulation. Figure [TO] presents snapshots of the vertical velocity at the solar surface for 
simulations with grid spacings of 100, 50, 25, and 12.5 km. The results show the granulation 
structure, which has a characteristic size of about 1 Mm and consists of relatively slow upflows 
occupying most of the area and fast downdrafts concentrated in the intergranular lanes. The 
primary effect of the decreasing grid spacing is in resolving the small-scale structures and 
dynamics in the intergranular lanes; the primary granulation structure is well-resolved at 
all resolutions. This is also evident from the turbulent energy spectra shown in Figure [TTJ 
The energy spectra show that resolving small-scale turbulence may significantly reduce the 
large-scale energy of the granulation. This effect has to be taken into account in stellar 
simulations which are generally performed with relatively low resolution. 

The small-scale dynamics of the intergranular lanes of solar convection, illustrated in 
Figure [l2l is associated with shearing flows and plays a critical role in the formation of 


tornado-like vortex tubes (Kitiashvili et al. 2012), excitation of acoustic waves (Kitiashvili 


et al. 2011), creation of local dynamos (Kitiashvili et al. 2013b), etc. For example, powerful 


vortex tubes are evident in Fig. [12] as rounded structures, ~ 100 km across, with dark cores 
corresponding to nearly supersonic downdrafts. Undoubtedly such shearing and twisting 
flows in the intergranular lanes are also extremely important on other stars. These effects 
will be a topic in future investigations using the StellarBox code. 


5.2. Structure of Granulation on the Sun and Moderate Mass Stars 

As stellar mass increases, the outer convection zone shrinks, and turbulent motions 
become more vigorous because of the increased energy flux. The granulation structure also 
changes quite significantly. Figure [13] shows the distribution of the vertical velocity on the 
surface of the Sun and five main-sequence stars with masses of 1.17 M 0 , 1.29 M 0 , 1.35 M 0 , 
1.47 M 0 , and 1.60 M 0 . The simulations were performed for a computational domain of 
100 x 100 Mm horizontally, 40 — 50 Mm in depth, and 1 Mm above the photosphere. The 
initial conditions in each case are for standard zero-age main-sequence models calculated 
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for the solar composition using the CESAM code (Morel 1997). The initial conditions for 


the solar simulations are constructed by combining the standard model S (([Christensen- 


Dalsgaard et al.||1996 )) and the VAL model of the solar atmosphere ((Vernazza et al.||1973[)). 


These hydrostatic initial conditions are perturbed by random velocity fluctuations to initiate 
convective motion, and simulation runs are continued for several hours of stellar time until 
statistically stationary conditions are achieved. No magnetic fields were included in these 
simulations. 


The simulation results show that the characteristic size of granulation increases from 
~ 1 Mm in the case of the Sun to more than 10 Mm for the 1.60 M 0 A-type star. For this 
latter type of star, standard stellar evolution theory using a mixing-length model predicts 
that the outer convection develops only in the hydrogen and helium ionization zones, sepa¬ 
rated by a stable layer. However, the simulations show that these layers are mixed, creating a 
shallow convection zone with a well-defined granulation pattern. For reference, the pressure 
scale heights for the Sun and the 1.17 M 0 , 1.29 M 0 , 1.35 M 0 , 1.47 M 0 , and 1.60 M 0 stars 
are approximately 140, 173, 236, 267, 270, and 359 km, respectively. 


An interesting question concerns the large-scale organization of granulation, observed 
on the Sun as meso- and super-granulation. Our results show that the solar granulation 
tends to form rather irregular clusters of several granules. However, the simulations do not 
explicitly show any large-scale (20-30 Mm) patterns which could correspond to supergran¬ 
ulation, at least for zero rotation and the simulation domain dimensions tried so far (the 
domain extended to 20 Mm deep for stars with masses < 1.35 M 0 and to approximately 
50 Mm deep for heavier stars). This means that some essential physics may be missing in 
these simulations, e.g., large-scale magnetic fields or rotation. This question requires further 
investigation. Mesoscale clustering can be seen in the more massive stars. It is particu¬ 
larly pronounced in the 1.47 M 0 case, in which large-scale convection cells become crossed, 
“shredded”, by a series of aligned intergranular lanes. This type of granular instability also 
exists in solar granulation (Kitiashvili et al. 2012) but is substantially less pronounced. The 
shredding process is accompanied by generation of intense acoustic waves which are seen as 
diffuse darker patches in the snapshots. 


These results demonstrate the capabilities of the StellarBox code in simulating turbulent 
stellar convection and reveal very complex multi-scale convective structures. A more detailed 
understanding of their dynamics will be a goal of our future studies. 
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5.3. Magnetic Field Structuring 


To illustrate the effects of magnetic fields on near-surface stellar convection, we present 
simulation results following the injection of a uniform 100 G vertical magnetic field into a fully 
developed convection layer. The boundary conditions conserve the mean vector magnetic 
field but do not prescribe any field structure. The domain is 12.8 x 12.8 Mm horizontally, 
5.5 Mm in depth below the photosphere, and 0.5 Mm above the photosphere. Results for 


the Sun and a 1.35 M 0 star are shown in Figure 14. In the case of the Sun, the initially 
uniform magnetic field becomes concentrated into compact, self-organized, 3-4 Mm wide 
pore-like magnetic structures maintained by strong dowdrafts. This process of spontaneous 


formation has been previously described in detail by Kitiashvili et al. (2010) for simulations 


in a smaller, 6x6 Mm domain. The new simulations confirm these results and show that 
the larger domain did not lead to formation of a larger structure. Instead, two compact 
structures of a similar size were formed, indicating that the structures are independent of 
the simulation domain size. What determines the scale of these self-organized magnetic 
domains has not been established. Curiously, the simulations for the 1.35 M 0 star did not 
result in formation of such large-scale structures, but instead the magnetic field became 
concentrated in small-scale patches in the intergranular lanes. In addition, the magnetic 
field formed diffuse patches of 100-200 G in the bodies of the stellar granules, whereas the 
magnetic fields inside the solar granules were much weaker. The reason for this difference 
is also unclear. Nevertheless, this indicates that the background (‘basal’) magnetic field 
may have quite different structures on different types of stars. From the numerical point of 
view, our simulations have shown that, for simulating the process of spontaneous structure 
formation, the grid spacing should be 25 km or smaller. For coarser grids the structures 
do not form. This indicates the importance of accurate simulation of the flow dynamics in 
intergranular lanes. 


6. Summary 

We have presented the basic features and some test results of the 3D radiative MHD 
code ‘StellarBox’. The code is designed to accurately simulate turbulent magnetoconvection 
processes in solar and stellar envelopes. It is based on a high-order Pade finite-difference 
scheme and implements subgrid-scale Large-Eddy Simulation (LES) turbulence modeling. 
This is expected to provide a more accurate description, as compared to using more ad hoc 
turbulence models or lower-accuracy numerics, of complex physical phenomena in the highly 
turbulent radiating plasma in the Sun and other stars. Such phenomena include the for¬ 
mation and dynamics of surface granulation, large-scale convective structures, excitation 
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of acoustic, gravity, and MHD waves, magnetic self-organization, etc. The code has been 
carefully tested using standard CFD and MHD solutions and shows good accuracy and 
robustness. 

The code’s capabilities were demonstrated by performing simulations of the upper con¬ 
vection zones of the Sun and several moderate mass stars, first without magnetic field and 
then with an imposed, initially uniform, vertical magnetic field. The solar convection simu¬ 
lations were performed for different numerical grid spacings, from 6.25 km to 100 km. The 
results show that, while the coarse, 100 km grid is sufficient for resolving the granulation 
structure, the rich dynamics of intergranular lanes can only be resolved with smaller grid 
spacings, e.g. 12 km. The importance of such resolution is revealed by the observation that 
the intergranular lanes are a source of strong, almost supersonic shearing flows and compact 
(~ 100 km across) vortex tubes, which play important roles in the formation of magnetic 
structures. 

For an initial comparative analysis of stellar convection, we performed large-scale (100 x 
100 x 40 Mm) simulations for six main-sequence stars with masses from 1.0 M 0 to 1.60 M 0 . 
The results showed a substantial increase in the granulation size with stellar mass, from 1 
to 20 Mm. Further, the granulation for stars more massive than the sun was often clustered 
on mesogranulation scales (on the scale of several granules). 

It was found that magnetic field effects can be quite different in stars of different classes. 
The simulations showed that, in the solar mass case, an initially uniform magnetic field forms 
self-organized, stable ‘pore-like’ structures of the size of several granules. But, in the case of 
more massive stars, such structures are not formed, and the magnetic field is distributed in 
the intergranular lanes in the form of small kilogauss magnetic flux tubes. 

In conclusion, 3D radiative MHD simulations, which are becoming more and more feasi¬ 
ble on modern supercomputer systems, allow us to simulate stellar magnetoconvection with a 
great degree of realism and provide an important tool for understanding the complex physics 
of turbulent stellar envelopes. 
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Fig. 1.— Illustration of the block scheme of the StellarBox code. 
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Fig. 2.— The 18 rays used to cover direction space O. 







(sec. / step / processor / grid point) 



Fig. 3.— Scaling results for StellarBox: a, time per step per processor per grid point, as a 
function of the number of processors. 
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Fig. 4.— Sod shock-tube test: a) Pressure and b) density profiles at t — 0.1. Solid line is 
the numerical solution and symbols represent the exact Riemann solver solution (ExactRS). 





Fig. 5.— Solution of the Orszag-Tang problem at t = 0.25, obtained with the StellarBox 
using the Pade scheme: a) density (the linear color map ranges from 0.06 to 0.52), b) 
magnetic energy (from 2.xl0 -8 to 0.3), c) kinetic energy (from 2.xl0 -8 to 0.65). 
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Fig. 6.— Solution of the Orszag-Tang problem at t — 0.25, obtained with StellarBox: a) 
density and b) pressure profiles y = 0.4277. 



Fig. 7.— Brio and Wu shock-tube test at t = 0.1: the solid line is the numerical result 
obtained with the Pade scheme, and the dashed line is the solution provided by the exact 
MHD Riemann solver of Torrilhon (2003): a) density profile, b) pressure profile, c) magnetic 
filed profile. 
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Fig. 8.— Comparison of radiative transfer calculations over the full domain; red line with 
squares: full ray-tracing; black line: diffusive radiation treatment in the optically thick region 
(below 2 : = —2.277 Mm). 
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Fig. 9.— Comparison of radiative transfer calculations below z = —2 Mm: same curves as 
in Fig. [8j 
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Fig. 10.— Snapshots of the simulated solar granulation (vertical velocity shown) at the 
photosphere for different resolutions: a) 100 km, b ) 50 km, c) 25 km and d) 12.5 km. 
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Fig. 13.— Variation of the scales of granulation for main-sequence stars with increasing 
stellar mass; from top, left to right: 1 M Q , 1.17 M 0 , 1.29 M Q , 1.35 M 0 , 1.47 M 0 , 1.60 M 0 . 
Distribution of the vertical velocity is plotted for a range of ±6 km/s in all plots. 
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Fig. 14.— Distribution of magnetic field in the photosphere from simulations with an initially 
uniform vertical 100 G magnetic held: a) in the Sun, and b) in a 1.35M 0 star. Panels c) 
and d) show corresponding vertical cuts along the x-axis at the locations indicated by the 
arrows in a) and b). The fields shown are typical of the statistically stationary state. 





















